import numpy as np
import matplotlib.pyplot as plt
# 关键修改：使用Windows默认中文支持字体（无需额外安装）
plt.rcParams["font.family"] = ["SimHei", "Microsoft YaHei", "FangSong", "KaiTi"]
plt.rcParams["axes.unicode_minus"] = False  # 解决负号显示问题

def generate_glv_communities(num_samples=10, N=100, t_end=100, dt=0.01,
                             r_range=(-0.1, 0.5), alpha_mag=0.01,
                             abundance_thresh=0.001):
    communities_s = []
    communities_p = []
    t = np.arange(0, t_end, dt)

    for _ in range(num_samples):
        # 优化参数范围，避免数值爆炸
        r = np.random.uniform(r_range[0], r_range[1], size=N)
        alpha = np.random.normal(0, alpha_mag, size=(N, N))
        diag = -np.abs(np.diag(alpha)) - 0.1  # 增强种内竞争，提升稳定性
        np.fill_diagonal(alpha, diag)
        x0 = np.random.uniform(0.01, 0.1, size=N)

        # 求解GLV方程，添加数值截断
        x = np.zeros((len(t), N))
        x[0] = x0
        for i in range(1, len(t)):
            dxdt = x[i - 1] * (r + np.dot(alpha, x[i - 1]))
            x[i] = x[i - 1] + dxdt * dt
            x[i] = np.clip(x[i], 0, 1e3)  # 限制丰度上限，防止溢出
            x[i] = np.maximum(x[i], 0)    # 丰度非负

        # 处理最终状态
        x_final = x[-1]
        s = (x_final > abundance_thresh).astype(int)
        total = np.sum(x_final)
        p = x_final / total if total > 0 else np.zeros(N)

        communities_s.append(s)
        communities_p.append(p)

    return np.array(communities_s), np.array(communities_p)

# 生成样本
num_samples = 10
N = 100
communities_s, communities_p = generate_glv_communities(
    num_samples, N,
    r_range=(-0.1, 0.3),
    alpha_mag=0.01
)

# 输出信息
print(f"生成了{num_samples}个群落样本，物种库大小为{N}")
print(f"物种存在矩阵形状: {communities_s.shape}")
print(f"物种丰度矩阵形状: {communities_p.shape}")

# 第一个样本详情
print("\n第一个样本信息:")
exist_count = np.sum(communities_s[0])
print("存在的物种数量:", exist_count)
if exist_count > 0:
    top5_idx = np.argsort(communities_p[0])[::-1][:5]
    print("丰度前5的物种及其丰度:")
    for idx in top5_idx:
        if communities_s[0][idx] == 1:
            print(f"物种{idx}: 丰度={communities_p[0][idx]:.4f}")
else:
    print("所有物种丰度均低于阈值，无存在物种")

# 可视化
plt.figure(figsize=(10, 4))
plt.bar(range(N), communities_p[0], alpha=0.7, color='#1f77b4')
plt.xlabel("物种索引")
plt.ylabel("丰度占比")
plt.title("第一个群落样本的物种丰度分布")
plt.grid(axis='y', alpha=0.3)
plt.show()